R Notebook

#Importing Data

#Creating a function to plot the average value of a specified variable in one year over time.

generalplotfunction=function(z){
  data=new_lakedata%>%
    filter(characteristic_name==z)%>%
    group_by(year)%>%
    summarise(mean=mean(result_value))
graph=ggplot(data,aes(x=year,y=mean))+
    geom_line()+ geom_point()+geom_vline(xintercept = 2002,color = "blue", linewidth=1.5)+geom_vline(xintercept = 2008,  color = "blue", size=1.5)+labs(title=z)
return(graph)
}

#Plot all the variables

#target <- c(unique(lakedata$characteristic_name))
#significant_variables=filter(lakedata, characteristic_name %in% target)
#average_sigvariables_by_year=significant_variables%>%
#  group_by(year,characteristic_name)%>%
 # summarise(Mean=mean(result_value))
#ggplot(average_sigvariables_by_year,aes(x=year,y=Mean))+
#  geom_line()+ geom_point()+geom_vline(xintercept = 2002,color = "blue", size=1.5)+geom_vline(xintercept = 2008,  color = "blue", size=1.5)+facet_wrap(vars(characteristic_name),scales = "free_y")

characteristic.names=c(unique(lakedata$characteristic_name))
plot_list2=list()
for( i in 1:length(characteristic.names)){
   plot.data= generalplotfunction(characteristic.names[i])
 plot_list2[[i]] <- plot.data
}

Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
ℹ Please use `linewidth` instead.
This warning is displayed once every 8 hours.
Call `lifecycle::last_lifecycle_warnings()` to see where this warning was generated.

plot_list2

[[1]]

[[2]]

[[3]]

[[4]]

[[5]]

[[6]]

[[7]]

[[8]]

[[9]]

[[10]]

[[11]]

[[12]]

[[13]]

[[14]]

[[15]]

[[16]]

[[17]]

[[18]]

[[19]]

[[20]]

[[21]]

[[22]]
`geom_line()`: Each group consists of only one observation.
ℹ Do you need to adjust the group aesthetic?

[[23]]

#Plot the six main variables

library(ggpubr)

PARTP_timeseries <- generalplotfunction("PARTP")
PARTN_timeseries <- generalplotfunction("PARTN")
TDP_timeseries <- generalplotfunction("TDP")
SO4_timeseries <- generalplotfunction("SO4")
TDN_timeseries <- generalplotfunction("TDN")
CHLA_timeseries <- generalplotfunction("CHLA")

ggarrange(PARTP_timeseries, PARTN_timeseries, TDP_timeseries,
          SO4_timeseries, TDN_timeseries, CHLA_timeseries+ rremove("x.text"), 
          ncol = 2, nrow = 3)

# confused by this part of the code (Karen)
target2=unique(new_lakedata$characteristic_name)
target2.vec=as.data.frame(target2)
 class(target2)

[1] "character"

 
for (i in target2) {
  generalplotfunction(i)
}

#Plotting Histograms

histogramfunction= function(charateristic)
{
  datafunc=new_lakedata%>%
    filter(characteristic_name==charateristic)%>%
    select(result_value)
  graph=ggplot(datafunc, aes(x = result_value)) +
    geom_histogram() +labs(title=charateristic)
  return(graph)
}

PARTP_hist <- histogramfunction("PARTP")
PARTN_hist <- histogramfunction("PARTN")
CHLA_hist <- histogramfunction("CHLA")
TDP_hist <- histogramfunction("TDP")
SO4_hist <- histogramfunction("SO4")
TDN_hist <- histogramfunction("TDN")

ggarrange(PARTP_hist, PARTN_hist, CHLA_hist, TDP_hist, SO4_hist,
          TDN_hist+ rremove("x.text"), 
          ncol = 2, nrow = 3)

`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

# confused by this code here (Karen)
histogramfunction(target2[1])

`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

for (i in target2) {
  histogramfunction(target2[i])
}

# creating a function that assesses the data's normality with shapiro test
normalfunction <- function(names) {
  n.data=new_lakedata %>%
    filter(characteristic_name==names) %>%
    select(result_value)
   num.data=as.numeric(n.data$result_value)
  norm.model=shapiro.test(num.data)
 return(norm.model$p.value)
}

# histograms for every variable 
characteristic.names=c(unique(new_lakedata$characteristic_name))
plot_list.norm <- list()
for( i in 1:length(characteristic.names)){
   plot.data <- histogramfunction(characteristic.names[i])
 plot_list.norm[[i]] <- plot.data
}
plot_list.norm

[[1]]
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

[[2]]
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

[[3]]
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

[[4]]
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

[[5]]
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

[[6]]
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

[[7]]
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

[[8]]
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

[[9]]
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

[[10]]
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

[[11]]

[[12]]
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

[[13]]
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

[[14]]
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

[[15]]
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

[[16]]
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

[[17]]
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

[[18]]
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

[[19]]
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

[[20]]
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

[[21]]
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

[[22]]
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

[[23]]
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

#Phosphorus models

#No need these just running mixed models below
library(lme4)

Loading required package: Matrix

Attaching package: ‘Matrix’

The following objects are masked from ‘package:tidyr’:

    expand, pack, unpack

library(lmerTest)

Attaching package: ‘lmerTest’

The following object is masked from ‘package:lme4’:

    lmer

The following object is masked from ‘package:stats’:

    step

library(MuMIn)

Partpdata <- lakedata_no_recovery %>%
  filter(characteristic_name=="PARTP")

m1 <- aov(result_value~Treatment+as.factor(year),data = Partpdata)
m2 <- glm(result_value~Treatment,data = Partpdata)
P_model_MM <- lmer(result_value~Treatment+(1|as.factor(Partpdata$year)), data=Partpdata, REML=FALSE)
summary(P_model_MM)

Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: result_value ~ Treatment + (1 | as.factor(Partpdata$year))
   Data: Partpdata

     AIC      BIC   logLik deviance df.resid 
   793.9    806.5   -392.9    785.9      170 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-2.2084 -0.3887 -0.0675  0.3421  6.2048 

Random effects:
 Groups                    Name        Variance Std.Dev.
 as.factor(Partpdata$year) (Intercept) 0.9751   0.9875  
 Residual                              4.8316   2.1981  
Number of obs: 174, groups:  as.factor(Partpdata$year), 17

Fixed effects:
                Estimate Std. Error      df t value Pr(>|t|)    
(Intercept)       2.9024     0.3750 22.0956   7.739 9.90e-08 ***
TreatmentDuring   3.8171     0.6339 17.4030   6.021 1.24e-05 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr)
TretmntDrng -0.592

ggplot(data=Partpdata, aes(x=Treatment,y=result_value))+
  geom_boxplot()

Tdpdata <- lakedata_no_recovery %>%
  filter(characteristic_name=="TDP")

summary(aov(result_value~Treatment,data = Tdpdata))

             Df Sum Sq Mean Sq F value   Pr(>F)    
Treatment     1   52.4   52.36   29.08 2.21e-07 ***
Residuals   176  316.9    1.80                     
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

summary(glm(result_value~Treatment,data = Tdpdata))

Call:
glm(formula = result_value ~ Treatment, data = Tdpdata)

Coefficients:
                Estimate Std. Error t value Pr(>|t|)    
(Intercept)       2.1667     0.1329  16.307  < 2e-16 ***
TreatmentDuring   1.0965     0.2033   5.393 2.21e-07 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for gaussian family taken to be 1.800588)

    Null deviance: 369.26  on 177  degrees of freedom
Residual deviance: 316.90  on 176  degrees of freedom
AIC: 613.81

Number of Fisher Scoring iterations: 2

P_model_MM <- lmer(result_value~Treatment+(1|as.factor(Tdpdata$year)), data=Tdpdata, REML=FALSE)
summary(P_model_MM)

Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: result_value ~ Treatment + (1 | as.factor(Tdpdata$year))
   Data: Tdpdata

     AIC      BIC   logLik deviance df.resid 
   598.5    611.3   -295.3    590.5      174 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-1.8745 -0.6243 -0.0480  0.4036  4.5137 

Random effects:
 Groups                  Name        Variance Std.Dev.
 as.factor(Tdpdata$year) (Intercept) 0.3387   0.582   
 Residual                            1.4442   1.202   
Number of obs: 178, groups:  as.factor(Tdpdata$year), 17

Fixed effects:
                Estimate Std. Error      df t value Pr(>|t|)    
(Intercept)       2.2200     0.2137 17.3158  10.389 7.33e-09 ***
TreatmentDuring   0.9889     0.3647 13.9181   2.711   0.0169 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr)
TretmntDrng -0.586

ggplot(data=Tdpdata, aes(x=Treatment,y=result_value))+
  geom_boxplot()

AICc(m1, m2,P_model_MM)

NA

#Nitrogen models

#No need these just running mixed models below
Partndata <- lakedata_no_recovery %>%
  filter(characteristic_name=="PARTN")

summary(aov(result_value~Treatment,data = Partndata))

             Df Sum Sq Mean Sq F value Pr(>F)    
Treatment     1  85667   85667   221.1 <2e-16 ***
Residuals   202  78263     387                   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

summary(glm(result_value~Treatment,data = Partndata))

Call:
glm(formula = result_value ~ Treatment, data = Partndata)

Coefficients:
                Estimate Std. Error t value Pr(>|t|)    
(Intercept)        40.56       1.74   23.32   <2e-16 ***
TreatmentDuring    42.38       2.85   14.87   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for gaussian family taken to be 387.442)

    Null deviance: 163931  on 203  degrees of freedom
Residual deviance:  78263  on 202  degrees of freedom
AIC: 1798.7

Number of Fisher Scoring iterations: 2

P_model_MM <- lmer(result_value~Treatment+(1|as.factor(Partndata$year)), data=Partndata, REML=FALSE)
summary(P_model_MM)

Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: result_value ~ Treatment + (1 | as.factor(Partndata$year))
   Data: Partndata

     AIC      BIC   logLik deviance df.resid 
  1798.6   1811.9   -895.3   1790.6      200 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-2.3521 -0.5778 -0.0526  0.3303  6.0707 

Random effects:
 Groups                    Name        Variance Std.Dev.
 as.factor(Partndata$year) (Intercept)  18.31    4.279  
 Residual                              365.69   19.123  
Number of obs: 204, groups:  as.factor(Partndata$year), 17

Fixed effects:
                Estimate Std. Error     df t value Pr(>|t|)    
(Intercept)       40.610      2.174 16.486   18.68 1.59e-12 ***
TreatmentDuring   42.057      3.645 13.784   11.54 1.82e-08 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr)
TretmntDrng -0.596

ggplot(data=Partndata, aes(x=Treatment,y=result_value))+
  geom_boxplot()

TdNdata <- lakedata_no_recovery %>%
  filter(characteristic_name=="TDN")

summary(aov(result_value~Treatment,data = TdNdata))

             Df Sum Sq Mean Sq F value   Pr(>F)    
Treatment     1  59161   59161   40.01 2.02e-09 ***
Residuals   176 260245    1479                     
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

summary(glm(result_value~Treatment,data = TdNdata))

Call:
glm(formula = result_value ~ Treatment, data = TdNdata)

Coefficients:
                Estimate Std. Error t value Pr(>|t|)    
(Intercept)      215.775      3.807  56.672  < 2e-16 ***
TreatmentDuring   36.857      5.827   6.325 2.02e-09 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for gaussian family taken to be 1478.668)

    Null deviance: 319406  on 177  degrees of freedom
Residual deviance: 260245  on 176  degrees of freedom
AIC: 1808.3

Number of Fisher Scoring iterations: 2

P_model_MM <- lmer(result_value~Treatment+(1|as.factor(TdNdata$year)), data=TdNdata, REML=FALSE)
summary(P_model_MM)

Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: result_value ~ Treatment + (1 | as.factor(TdNdata$year))
   Data: TdNdata

     AIC      BIC   logLik deviance df.resid 
  1790.0   1802.7   -891.0   1782.0      174 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-3.8432 -0.5816 -0.0552  0.4335  4.5410 

Random effects:
 Groups                  Name        Variance Std.Dev.
 as.factor(TdNdata$year) (Intercept)  268.1   16.37   
 Residual                            1167.1   34.16   
Number of obs: 178, groups:  as.factor(TdNdata$year), 17

Fixed effects:
                Estimate Std. Error      df t value Pr(>|t|)    
(Intercept)      215.768      6.033  20.436   35.76  < 2e-16 ***
TreatmentDuring   36.013     10.290  16.456    3.50  0.00286 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr)
TretmntDrng -0.586

ggplot(data=TdNdata, aes(x=Treatment,y=result_value))+
  geom_boxplot()

NA
NA

#Permutation

#Inital Permutation code
permutation_function <- function(){
  reshuffled <- Partpdata
  reshuffled$TreatmentNEW <- sample(reshuffled$Treatment, size = nrow(reshuffled), replace = F)
  
  before <- reshuffled %>% subset(TreatmentNEW == "Before") %>% summarise(mean = mean(result_value))
  during <- reshuffled %>% subset(TreatmentNEW == "During") %>% summarise(mean = mean(result_value))

  return(before-during)
}

gaps <- c()

for (i in 1:1000){
  gaps[i] <- permutation_function()
}

hist(unlist(gaps))

#Creating function for LMM

  mixed_model_data <- lakedata_no_recovery %>%
    filter(characteristic_name == characteristic) %>%
    select(result_value, Treatment, year)

Error in `filter()`:
ℹ In argument: `characteristic_name == characteristic`.
Caused by error:
! object 'characteristic' not found
Backtrace:
 1. ... %>% select(result_value, Treatment, year)
 4. dplyr:::filter.data.frame(., characteristic_name == characteristic)
 5. dplyr:::filter_rows(.data, dots, by)
 6. dplyr:::filter_eval(...)
 8. mask$eval_all_filter(dots, env_filter)
 9. dplyr (local) eval()

#creating function for Permutation

library(tidyverse)
#Making Permutation Function
filterfunction=function(characteristic){
  data.perm=lakedata_no_recovery%>%
  filter(characteristic_name==characteristic)%>%
  select(Treatment,result_value,year)
  return(data.perm)
  }
permutation_function <- function(data12){
  reshuffled <- data12
  reshuffled$TreatmentNEW <- sample(reshuffled$Treatment, size = nrow(reshuffled), replace = F)
  model=glm(result_value~TreatmentNEW, data = reshuffled)
   #lmer(result_value~TreatmentNEW+(1|as.factor(reshuffled$year)), data=reshuffled, REML=FALSE)
  summary.model=summary(model)
  return(summary.model$coefficients[2,1])
}


final.perm.function= function(characteristic){
filtered.data=filterfunction(characteristic)
estimates <- c()
for (i in 1:10000){
  estimates[i] <- permutation_function(filtered.data)
}
estimates=as.data.frame(estimates)
 data.perm=lakedata_no_recovery%>%
  filter(characteristic_name==characteristic)%>%
  select(Treatment,result_value,year)
model.new=summary(glm(result_value~Treatment, data = data.perm))
lol=model.new$coefficients[2,1]
plot_s=ggplot(estimates, aes(x=estimates)) + 
 geom_histogram( colour="black", fill="white")+ geom_vline(aes(xintercept=lol),
            color="blue", linetype="dashed", size=1)+labs(title =characteristic)
return(plot_s)
}
characteristic.names=c(unique(lakedata_no_recovery$characteristic_name))
plot_list <- list()
for( i in 1:(length(characteristic.names)-1)){
   plot.data= final.perm.function(characteristic.names[i])
 plot_list[[i]] <- plot.data
}

plot_list

[[1]]
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

[[2]]
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

[[3]]
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

[[4]]
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

[[5]]
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

[[6]]
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

[[7]]
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

[[8]]
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

[[9]]
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

[[10]]
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

[[11]]
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

[[12]]
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

[[13]]
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

[[14]]
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

[[15]]
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

[[16]]
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

[[17]]
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

[[18]]
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

[[19]]
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

#Function for ALL LMMs

  pvalues[i]=  returnpvalueforllm(characteristic.names[i])

Error: couldn't evaluate grouping factor as.factor(datafortest$year) within model frame:error =Error in eval(substitute(makeFac(fac), list(fac = ff0)), frloc) : 
  object 'datafortest' not found
 Try adding grouping factor to data frame explicitly if possible

#Running Aarimas

library(forecast)

Registered S3 method overwritten by 'quantmod':
  method            from
  as.zoo.data.frame zoo 

Attaching package: ‘forecast’

The following object is masked from ‘package:ggpubr’:

    gghistogram

alldataforecastfunction=function(y){
  data1=lakedata%>%
    filter(characteristic_name==y)%>%
    select(result_value)
  model3=auto.arima(data1) 
  summary(model3)
  forecast_data3 <- forecast(model3, 20) 
  plot=plot(forecast_data3, main = y)
  return(plot)
}
beforeforecastfunction=function(y){
  data1=lakedata%>%
    filter(characteristic_name==y)%>%
    filter(year<2003)%>%
    select(result_value)
  model3=auto.arima(data1) 
  summary(model3)
  forecast_data3 <- forecast(model3, 20) 
  plot=plot(forecast_data3, main = y)
  return(plot)
}

duringforecastfunction=function(y){
  data1=lakedata%>%
    filter(characteristic_name==y)%>%
    filter(year<2007)%>%
    select(result_value)
  model3=auto.arima(data1) 
  summary(model3)
  forecast_data3 <- forecast(model3, 20) 
  plot=plot(forecast_data3, main = y)
  return(plot)
}
ariimalist=list()
for( i in 1:(length(characteristic.names)-1)){
   plot.data= alldataforecastfunction(characteristic.names[i])
 ariimalist[[i]] <- plot.data
}

beforeariimalist=list()
for( i in 1:(length(characteristic.names)-1)){
   plot.data= beforeforecastfunction(characteristic.names[i])
 beforeariimalist[[i]] <- plot.data
}

duringariimalist=list()
for( i in 1:(length(characteristic.names)-1)){
   plot.data= beforeforecastfunction(characteristic.names[i])
 duringariimalist[[i]] <- plot.data
}

ariimalist

[[1]]
[[1]]$mean
Time Series:
Start = 200 
End = 219 
Frequency = 1 
 [1] 389.2917 378.9795 372.7945 367.3273 364.4093 361.4367 360.1407 358.4647 357.9622 356.9704 356.8458 356.2234 356.2698 355.8536 355.9694 355.6737
[17] 355.8101 355.5889 355.7234 355.5513

[[1]]$lower
Time Series:
Start = 200 
End = 219 
Frequency = 1 
         80%      95%
200 368.5536 357.5755
201 351.9698 337.6717
202 343.1750 327.4953
203 335.7412 319.0206
204 331.8620 314.6325
205 327.9869 310.2797
206 326.2076 308.2446
207 323.9997 305.7550
208 323.1847 304.7746
209 321.8140 303.2034
210 321.4460 302.7065
211 320.5164 301.6143
212 320.3494 301.3344
213 319.6633 300.5052
214 319.5791 300.3152
215 319.0346 299.6390
216 318.9769 299.4786
217 318.5203 298.8973
218 318.4634 298.7391
219 318.0650 298.2209

[[1]]$upper
Time Series:
Start = 200 
End = 219 
Frequency = 1 
         80%      95%
200 410.0299 421.0080
201 405.9892 420.2872
202 402.4141 418.0937
203 398.9134 415.6341
204 396.9566 414.1861
205 394.8864 412.5936
206 394.0737 412.0368
207 392.9297 411.1743
208 392.7396 411.1497
209 392.1267 410.7374
210 392.2455 410.9850
211 391.9304 410.8325
212 392.1901 411.2052
213 392.0440 411.2020
214 392.3598 411.6237
215 392.3129 411.7085
216 392.6433 412.1416
217 392.6576 412.2805
218 392.9835 412.7078
219 393.0377 412.8817


[[2]]
[[2]]$mean
Time Series:
Start = 186 
End = 205 
Frequency = 1 
 [1] 4.778458 4.788540 4.794046 4.797054 4.798697 4.799594 4.800084 4.800352 4.800498 4.800578 4.800621 4.800645 4.800658 4.800665 4.800669 4.800671
[17] 4.800673 4.800673 4.800674 4.800674

[[2]]$lower
Time Series:
Start = 186 
End = 205 
Frequency = 1 
         80%      95%
186 4.313410 4.067229
187 4.222850 3.923392
188 4.177981 3.851856
189 4.148923 3.805823
190 4.126556 3.770745
191 4.107420 3.741006
192 4.090024 3.714141
193 4.073662 3.688976
194 4.057979 3.664914
195 4.042788 3.641638
196 4.027982 3.618972
197 4.013502 3.596813
198 3.999308 3.575098
199 3.985374 3.553784
200 3.971681 3.532841
201 3.958214 3.512245
202 3.944962 3.491976
203 3.931912 3.472018
204 3.919057 3.452357
205 3.906386 3.432978

[[2]]$upper
Time Series:
Start = 186 
End = 205 
Frequency = 1 
         80%      95%
186 5.243506 5.489687
187 5.354229 5.653687
188 5.410112 5.736237
189 5.445185 5.788285
190 5.470838 5.826648
191 5.491768 5.858182
192 5.510144 5.886027
193 5.527042 5.911728
194 5.543017 5.936083
195 5.558368 5.959518
196 5.573261 5.982271
197 5.587789 6.004477
198 5.602009 6.026218
199 5.615957 6.047546
200 5.629658 6.068498
201 5.643128 6.089098
202 5.656383 6.109369
203 5.669434 6.129328
204 5.682291 6.148991
205 5.694962 6.168369


[[3]]
[[3]]$mean
Time Series:
Start = 258 
End = 277 
Frequency = 1 
 [1] 1.759686 1.789013 1.837400 1.837400 1.837400 1.837400 1.837400 1.837400 1.837400 1.837400 1.837400 1.837400 1.837400 1.837400 1.837400 1.837400
[17] 1.837400 1.837400 1.837400 1.837400

[[3]]$lower
Time Series:
Start = 258 
End = 277 
Frequency = 1 
           80%       95%
258 -0.7035347 -2.007485
259 -0.7355183 -2.071925
260 -0.7692404 -2.149113
261 -0.7802297 -2.165919
262 -0.7911730 -2.182656
263 -0.8020710 -2.199323
264 -0.8129241 -2.215921
265 -0.8237330 -2.232452
266 -0.8344981 -2.248916
267 -0.8452201 -2.265314
268 -0.8558993 -2.281646
269 -0.8665364 -2.297914
270 -0.8771318 -2.314118
271 -0.8876860 -2.330260
272 -0.8981995 -2.346339
273 -0.9086728 -2.362356
274 -0.9191062 -2.378313
275 -0.9295003 -2.394209
276 -0.9398555 -2.410046
277 -0.9501723 -2.425824

[[3]]$upper
Time Series:
Start = 258 
End = 277 
Frequency = 1 
         80%      95%
258 4.222906 5.526856
259 4.313544 5.649951
260 4.444041 5.823913
261 4.455030 5.840720
262 4.465974 5.857456
263 4.476872 5.874123
264 4.487725 5.890722
265 4.498534 5.907253
266 4.509299 5.923716
267 4.520021 5.940114
268 4.530700 5.956447
269 4.541337 5.972715
270 4.551933 5.988919
271 4.562487 6.005060
272 4.573000 6.021139
273 4.583473 6.037157
274 4.593907 6.053113
275 4.604301 6.069010
276 4.614656 6.084847
277 4.624973 6.100625


[[4]]
[[4]]$mean
Time Series:
Start = 192 
End = 211 
Frequency = 1 
 [1] 0.2253348 0.2268350 0.2257525 0.2251003 0.2242551 0.2234396 0.2226097 0.2217818 0.2209528 0.2201240 0.2192951 0.2184662 0.2176373 0.2168084
[15] 0.2159794 0.2151505 0.2143216 0.2134927 0.2126638 0.2118349

[[4]]$lower
Time Series:
Start = 192 
End = 211 
Frequency = 1 
          80%        95%
192 0.1192894 0.06315242
193 0.1203658 0.06400449
194 0.1181074 0.06112358
195 0.1171231 0.05996338
196 0.1159186 0.05856879
197 0.1147899 0.05727428
198 0.1136464 0.05596479
199 0.1125092 0.05466381
200 0.1113718 0.05336303
201 0.1102355 0.05206411
202 0.1091001 0.05076642
203 0.1079656 0.04947007
204 0.1068319 0.04817501
205 0.1056990 0.04688122
206 0.1045669 0.04558871
207 0.1034357 0.04429745
208 0.1023053 0.04300744
209 0.1011757 0.04171867
210 0.1000469 0.04043112
211 0.0989189 0.03914480

[[4]]$upper
Time Series:
Start = 192 
End = 211 
Frequency = 1 
          80%       95%
192 0.3313801 0.3875172
193 0.3333041 0.3896655
194 0.3333976 0.3903815
195 0.3330775 0.3902372
196 0.3325915 0.3899413
197 0.3320893 0.3896049
198 0.3315729 0.3892546
199 0.3310544 0.3888998
200 0.3305339 0.3885426
201 0.3300125 0.3881839
202 0.3294901 0.3878238
203 0.3289668 0.3874623
204 0.3284427 0.3870995
205 0.3279177 0.3867355
206 0.3273920 0.3863702
207 0.3268654 0.3860036
208 0.3263379 0.3856358
209 0.3258097 0.3852667
210 0.3252807 0.3848965
211 0.3247509 0.3845250


[[5]]
[[5]]$mean
Time Series:
Start = 204 
End = 223 
Frequency = 1 
 [1] 47.9064 47.9064 47.9064 47.9064 47.9064 47.9064 47.9064 47.9064 47.9064 47.9064 47.9064 47.9064 47.9064 47.9064 47.9064 47.9064 47.9064 47.9064
[19] 47.9064 47.9064

[[5]]$lower
Time Series:
Start = 204 
End = 223 
Frequency = 1 
         80%       95%
204 9.206788 -11.27955
205 9.206788 -11.27955
206 9.206788 -11.27955
207 9.206788 -11.27955
208 9.206788 -11.27955
209 9.206788 -11.27955
210 9.206788 -11.27955
211 9.206788 -11.27955
212 9.206788 -11.27955
213 9.206788 -11.27955
214 9.206788 -11.27955
215 9.206788 -11.27955
216 9.206788 -11.27955
217 9.206788 -11.27955
218 9.206788 -11.27955
219 9.206788 -11.27955
220 9.206788 -11.27955
221 9.206788 -11.27955
222 9.206788 -11.27955
223 9.206788 -11.27955

[[5]]$upper
Time Series:
Start = 204 
End = 223 
Frequency = 1 
         80%      95%
204 86.60602 107.0924
205 86.60602 107.0924
206 86.60602 107.0924
207 86.60602 107.0924
208 86.60602 107.0924
209 86.60602 107.0924
210 86.60602 107.0924
211 86.60602 107.0924
212 86.60602 107.0924
213 86.60602 107.0924
214 86.60602 107.0924
215 86.60602 107.0924
216 86.60602 107.0924
217 86.60602 107.0924
218 86.60602 107.0924
219 86.60602 107.0924
220 86.60602 107.0924
221 86.60602 107.0924
222 86.60602 107.0924
223 86.60602 107.0924


[[6]]
[[6]]$mean
Time Series:
Start = 403 
End = 422 
Frequency = 1 
 [1] 353.7947 362.8866 345.6985 362.1347 346.4175 361.4472 347.0749 360.8185 347.6762 360.2436 348.2259 359.7178 348.7286 359.2371 349.1883 358.7975
[17] 349.6087 358.3956 349.9931 358.0280

[[6]]$lower
Time Series:
Start = 403 
End = 422 
Frequency = 1 
         80%      95%
403 294.6680 263.3682
404 299.8022 266.4074
405 281.3067 247.2198
406 297.6157 263.4613
407 280.6945 245.9028
408 295.5756 260.7053
409 280.0890 244.6288
410 293.6636 258.1139
411 279.4849 243.3866
412 291.8640 255.6660
413 278.8782 242.1677
414 290.1638 253.3441
415 278.2664 240.9659
416 288.5517 251.1331
417 277.6477 239.7764
418 287.0181 249.0204
419 277.0210 238.5954
420 285.5548 246.9953
421 276.3858 237.4205
422 284.1546 245.0485

[[6]]$upper
Time Series:
Start = 403 
End = 422 
Frequency = 1 
         80%      95%
403 412.9215 444.2213
404 425.9711 459.3659
405 410.0902 444.1771
406 426.6538 460.8082
407 412.1405 446.9321
408 427.3188 462.1891
409 414.0609 449.5211
410 427.9734 463.5230
411 415.8675 451.9658
412 428.6231 464.8211
413 417.5736 454.2841
414 429.2719 466.0916
415 419.1909 456.4913
416 429.9226 467.3412
417 420.7290 458.6003
418 430.5769 468.5747
419 422.1964 460.6220
420 431.2363 469.7959
421 423.6003 462.5656
422 431.9013 471.0075


[[7]]
[[7]]$mean
Time Series:
Start = 217 
End = 236 
Frequency = 1 
 [1] 475.6166 452.2344 448.2061 472.6186 480.7788 483.5064 484.4181 484.7229 484.8248 484.8588 484.8702 484.8740 484.8753 484.8757 484.8758 484.8759
[17] 484.8759 484.8759 484.8759 484.8759

[[7]]$lower
Time Series:
Start = 217 
End = 236 
Frequency = 1 
         80%      95%
217 384.3057 335.9686
218 347.5262 292.0971
219 337.3991 278.7414
220 359.3426 299.3778
221 367.2302 307.1212
222 369.9274 309.8023
223 370.8357 310.7088
224 371.1401 311.0130
225 371.2419 311.1148
226 371.2760 311.1488
227 371.2873 311.1602
228 371.2911 311.1640
229 371.2924 311.1653
230 371.2928 311.1657
231 371.2930 311.1659
232 371.2930 311.1659
233 371.2930 311.1659
234 371.2931 311.1659
235 371.2931 311.1659
236 371.2931 311.1659

[[7]]$upper
Time Series:
Start = 217 
End = 236 
Frequency = 1 
         80%      95%
217 566.9274 615.2645
218 556.9425 612.3717
219 559.0131 617.6707
220 585.8947 645.8595
221 594.3274 654.4364
222 597.0854 657.2105
223 598.0006 658.1275
224 598.3057 658.4328
225 598.4076 658.5347
226 598.4417 658.5688
227 598.4530 658.5802
228 598.4569 658.5840
229 598.4581 658.5853
230 598.4585 658.5857
231 598.4587 658.5858
232 598.4587 658.5859
233 598.4588 658.5859
234 598.4588 658.5859
235 598.4588 658.5859
236 598.4588 658.5859


[[8]]
[[8]]$mean
Time Series:
Start = 163 
End = 182 
Frequency = 1 
 [1] 0.4826691 0.5106029 0.4880757 0.5121071 0.4867426 0.5083343 0.4872397 0.5056626 0.4896151 0.5037094 0.4919697 0.5018713 0.4937423 0.5002978
[15] 0.4950476 0.4991367 0.4959959 0.4983487 0.4966415 0.4978387

[[8]]$lower
Time Series:
Start = 163 
End = 182 
Frequency = 1 
          80%       95%
163 0.4230685 0.3915179
164 0.4458813 0.4116197
165 0.4170092 0.3793889
166 0.4402899 0.4022721
167 0.4143897 0.3760884
168 0.4353415 0.3967014
169 0.4139315 0.3751246
170 0.4309779 0.3914422
171 0.4144893 0.3747201
172 0.4269324 0.3862892
173 0.4147398 0.3738568
174 0.4231504 0.3814781
175 0.4145596 0.3726428
176 0.4198151 0.3772102
177 0.4140475 0.3711687
178 0.4169810 0.3734903
179 0.4132527 0.3694511
180 0.4145675 0.3702164
181 0.4122147 0.3675218
182 0.4124694 0.3672777

[[8]]$upper
Time Series:
Start = 163 
End = 182 
Frequency = 1 
          80%       95%
163 0.5422697 0.5738204
164 0.5753246 0.6095861
165 0.5591421 0.5967624
166 0.5839244 0.6219422
167 0.5590954 0.5973967
168 0.5813272 0.6199673
169 0.5605478 0.5993548
170 0.5803474 0.6198831
171 0.5647409 0.6045101
172 0.5804863 0.6211295
173 0.5691996 0.6100827
174 0.5805922 0.6222645
175 0.5729249 0.6148417
176 0.5807805 0.6233854
177 0.5760477 0.6189265
178 0.5812925 0.6247831
179 0.5787391 0.6225406
180 0.5821299 0.6264810
181 0.5810684 0.6257613
182 0.5832080 0.6283998


[[9]]
[[9]]$mean
Time Series:
Start = 198 
End = 217 
Frequency = 1 
 [1] 1.158958 1.149087 1.154031 1.154147 1.155066 1.155432 1.155733 1.155911 1.156031 1.156108 1.156158 1.156191 1.156212 1.156226 1.156235 1.156240
[17] 1.156244 1.156247 1.156248 1.156249

[[9]]$lower
Time Series:
Start = 198 
End = 217 
Frequency = 1 
          80%       95%
198 1.0733284 1.0279990
199 1.0515187 0.9998694
200 1.0456114 0.9882173
201 1.0383759 0.9770902
202 1.0331643 0.9686334
203 1.0283051 0.9610080
204 1.0239141 0.9541334
205 1.0197718 0.9477040
206 1.0158317 0.9416145
207 1.0120426 0.9357788
208 1.0083779 0.9301476
209 1.0048182 0.9246863
210 1.0013504 0.9193716
211 0.9979647 0.9141863
212 0.9946537 0.9091179
213 0.9914116 0.9041564
214 0.9882334 0.8992939
215 0.9851152 0.8945237
216 0.9820533 0.8898401
217 0.9790447 0.8852382

[[9]]$upper
Time Series:
Start = 198 
End = 217 
Frequency = 1 
         80%      95%
198 1.244587 1.289917
199 1.246654 1.298304
200 1.262451 1.319846
201 1.269919 1.331204
202 1.276968 1.341498
203 1.282560 1.349857
204 1.287552 1.357333
205 1.292050 1.364118
206 1.296231 1.370448
207 1.300174 1.376438
208 1.303939 1.382169
209 1.307564 1.387696
210 1.311074 1.393053
211 1.314487 1.398265
212 1.317816 1.403351
213 1.321069 1.408324
214 1.324255 1.413194
215 1.327378 1.417970
216 1.330443 1.422656
217 1.333454 1.427260


[[10]]
[[10]]$mean
Time Series:
Start = 161 
End = 180 
Frequency = 1 
 [1] 0.01977020 0.01517453 0.01515076 0.01515064 0.01515064 0.01515064 0.01515064 0.01515064 0.01515064 0.01515064 0.01515064 0.01515064 0.01515064
[14] 0.01515064 0.01515064 0.01515064 0.01515064 0.01515064 0.01515064 0.01515064

[[10]]$lower
Time Series:
Start = 161 
End = 180 
Frequency = 1 
            80%         95%
161 -0.01319113 -0.03063981
162 -0.02000168 -0.03862284
163 -0.02006018 -0.03869972
164 -0.02009258 -0.03874921
165 -0.02012482 -0.03879852
166 -0.02015703 -0.03884778
167 -0.02018921 -0.03889699
168 -0.02022136 -0.03894616
169 -0.02025348 -0.03899528
170 -0.02028557 -0.03904436
171 -0.02031763 -0.03909340
172 -0.02034966 -0.03914239
173 -0.02038167 -0.03919134
174 -0.02041364 -0.03924024
175 -0.02044559 -0.03928910
176 -0.02047751 -0.03933791
177 -0.02050940 -0.03938668
178 -0.02054126 -0.03943541
179 -0.02057309 -0.03948410
180 -0.02060490 -0.03953274

[[10]]$upper
Time Series:
Start = 161 
End = 180 
Frequency = 1 
           80%        95%
161 0.05273153 0.07018021
162 0.05035074 0.06897190
163 0.05036170 0.06900124
164 0.05039386 0.06905049
165 0.05042609 0.06909979
166 0.05045830 0.06914905
167 0.05049048 0.06919826
168 0.05052263 0.06924743
169 0.05055475 0.06929656
170 0.05058684 0.06934564
171 0.05061890 0.06939467
172 0.05065094 0.06944366
173 0.05068294 0.06949261
174 0.05071492 0.06954151
175 0.05074687 0.06959037
176 0.05077878 0.06963919
177 0.05081067 0.06968796
178 0.05084254 0.06973668
179 0.05087437 0.06978537
180 0.05090617 0.06983401


[[11]]
[[11]]$mean
Time Series:
Start = 217 
End = 236 
Frequency = 1 
 [1] 14.79651 14.60336 13.98920 13.98920 13.98920 13.98920 13.98920 13.98920 13.98920 13.98920 13.98920 13.98920 13.98920 13.98920 13.98920 13.98920
[17] 13.98920 13.98920 13.98920 13.98920

[[11]]$lower
Time Series:
Start = 217 
End = 236 
Frequency = 1 
           80%       95%
217  -6.262042 -17.40976
218 -10.546189 -23.85956
219 -11.574154 -25.10658
220 -11.608414 -25.15897
221 -11.642629 -25.21130
222 -11.676797 -25.26356
223 -11.710921 -25.31574
224 -11.744999 -25.36786
225 -11.779032 -25.41991
226 -11.813020 -25.47189
227 -11.846963 -25.52380
228 -11.880862 -25.57565
229 -11.914716 -25.62742
230 -11.948527 -25.67913
231 -11.982293 -25.73077
232 -12.016015 -25.78235
233 -12.049694 -25.83385
234 -12.083329 -25.88529
235 -12.116921 -25.93667
236 -12.150470 -25.98798

[[11]]$upper
Time Series:
Start = 217 
End = 236 
Frequency = 1 
         80%      95%
217 35.85506 47.00278
218 39.75291 53.06628
219 39.55256 53.08499
220 39.58682 53.13738
221 39.62104 53.18971
222 39.65521 53.24197
223 39.68933 53.29415
224 39.72341 53.34627
225 39.75744 53.39832
226 39.79143 53.45030
227 39.82537 53.50221
228 39.85927 53.55406
229 39.89312 53.60583
230 39.92693 53.65754
231 39.96070 53.70918
232 39.99442 53.76076
233 40.02810 53.81226
234 40.06174 53.86370
235 40.09533 53.91508
236 40.12888 53.96639


[[12]]
[[12]]$mean
Time Series:
Start = 254 
End = 273 
Frequency = 1 
 [1] 455.5767 476.4266 483.1584 485.3319 486.0336 486.2602 486.3334 486.3570 486.3646 486.3671 486.3679 486.3681 486.3682 486.3682 486.3683 486.3683
[17] 486.3683 486.3683 486.3683 486.3683

[[12]]$lower
Time Series:
Start = 254 
End = 273 
Frequency = 1 
         80%       95%
254 249.7012 140.71733
255 253.4953 135.48268
256 255.3964 134.82647
257 255.0293 133.11446
258 253.7772 130.82802
259 252.2305 128.34261
260 250.5948 125.80235
261 248.9381 123.25619
262 247.2825 120.72004
263 245.6348 118.19883
264 243.9972 115.69396
265 242.3703 113.20563
266 240.7540 110.73369
267 239.1483 108.27789
268 237.5529 105.83793
269 235.9676 103.41352
270 234.3923 101.00435
271 232.8269  98.61015
272 231.2710  96.23065
273 229.7245  93.86557

[[12]]$upper
Time Series:
Start = 254 
End = 273 
Frequency = 1 
         80%      95%
254 661.4521 770.4360
255 699.3578 817.3704
256 710.9204 831.4903
257 715.6345 837.5493
258 718.2901 841.2393
259 720.2900 844.1778
260 722.0719 846.8644
261 723.7758 849.4578
262 725.4468 852.0092
263 727.0994 854.5353
264 728.7385 857.0418
265 730.3660 859.5306
266 731.9824 862.0027
267 733.5882 864.4586
268 735.1837 866.8986
269 736.7689 869.3230
270 738.3442 871.7322
271 739.9097 874.1264
272 741.4655 876.5059
273 743.0120 878.8709


[[13]]
[[13]]$mean
Time Series:
Start = 250 
End = 269 
Frequency = 1 
 [1] 50.53113 49.67078 49.37097 49.26649 49.23008 49.21740 49.21297 49.21143 49.21090 49.21071 49.21064 49.21062 49.21061 49.21061 49.21061 49.21061
[17] 49.21061 49.21061 49.21061 49.21061

[[13]]$lower
Time Series:
Start = 250 
End = 269 
Frequency = 1 
         80%       95%
250 27.15526 14.780816
251 23.99225 10.398855
252 22.95790  8.975664
253 22.44936  8.253232
254 22.10161  7.740671
255 21.80857  7.299214
256 21.53581  6.884400
257 21.27163  6.481196
258 21.01199  6.084384
259 20.75542  5.692098
260 20.50139  5.303632
261 20.24968  4.918681
262 20.00016  4.537086
263 19.75277  4.158735
264 19.50745  3.783540
265 19.26413  3.411421
266 19.02278  3.042301
267 18.78334  2.676109
268 18.54576  2.312776
269 18.31002  1.952237

[[13]]$upper
Time Series:
Start = 250 
End = 269 
Frequency = 1 
         80%      95%
250 73.90700 86.28144
251 75.34931 88.94271
252 75.78404 89.76627
253 76.08362 90.27975
254 76.35855 90.71950
255 76.62622 91.13558
256 76.89014 91.54155
257 77.15123 91.94167
258 77.40981 92.33741
259 77.66600 92.72932
260 77.91990 93.11766
261 78.17156 93.50256
262 78.42106 93.88414
263 78.66845 94.26249
264 78.91377 94.63768
265 79.15709 95.00980
266 79.39844 95.37892
267 79.63788 95.74511
268 79.87545 96.10844
269 80.11120 96.46898


[[14]]
[[14]]$mean
Time Series:
Start = 216 
End = 235 
Frequency = 1 
 [1] 2.894592 3.350959 3.583770 3.702536 3.763123 3.794031 3.809799 3.817842 3.821945 3.824039 3.825107 3.825651 3.825929 3.826071 3.826143 3.826180
[17] 3.826199 3.826209 3.826214 3.826216

[[14]]$lower
Time Series:
Start = 216 
End = 235 
Frequency = 1 
             80%       95%
216  0.221504928 -1.193542
217  0.245602139 -1.398274
218  0.314228624 -1.416562
219  0.350581446 -1.423836
220  0.359004072 -1.443028
221  0.350476255 -1.472432
222  0.332687000 -1.507985
223  0.310086571 -1.546807
224  0.285082520 -1.587220
225  0.258938987 -1.628311
226  0.232309677 -1.669602
227  0.205529334 -1.710848
228  0.178767776 -1.751923
229  0.152110018 -1.792767
230  0.125597611 -1.833353
231  0.099249897 -1.873668
232  0.073074905 -1.913709
233  0.047074934 -1.953478
234  0.021249414 -1.992977
235 -0.004403641 -2.032211

[[14]]$upper
Time Series:
Start = 216 
End = 235 
Frequency = 1 
         80%      95%
216 5.567678 6.982725
217 6.456315 8.100191
218 6.853310 8.584101
219 7.054490 8.828907
220 7.167242 8.969274
221 7.237586 9.060494
222 7.286910 9.127582
223 7.325598 9.182491
224 7.358808 9.231111
225 7.389139 9.276388
226 7.417904 9.319815
227 7.445773 9.362150
228 7.473091 9.403782
229 7.500032 9.444910
230 7.526689 9.485640
231 7.553111 9.526028
232 7.579323 9.566107
233 7.605342 9.605895
234 7.631178 9.645404
235 7.656836 9.684643


[[15]]
[[15]]$mean
Time Series:
Start = 253 
End = 272 
Frequency = 1 
 [1] 6.680631 6.717498 6.986455 6.913629 7.080272 6.984759 7.105291 7.013643 7.108563 7.027735 7.105453 7.036224 7.100903 7.042256 7.096438 7.046964
[17] 7.092472 7.050805 7.089065 7.053996

[[15]]$lower
Time Series:
Start = 253 
End = 272 
Frequency = 1 
         80%      95%
253 6.373580 6.211037
254 6.382997 6.205924
255 6.622486 6.429813
256 6.546516 6.352178
257 6.703134 6.503490
258 6.606928 6.406917
259 6.721927 6.518987
260 6.629901 6.426760
261 6.720815 6.515554
262 6.639608 6.434145
263 6.714035 6.506831
264 6.644349 6.436904
265 6.706168 6.497208
266 6.646965 6.437711
267 6.698588 6.487979
268 6.648460 6.437504
269 6.691638 6.479449
270 6.649224 6.436641
271 6.685337 6.471617
272 6.649440 6.435280

[[15]]$upper
Time Series:
Start = 253 
End = 272 
Frequency = 1 
         80%      95%
253 6.987682 7.150225
254 7.051998 7.229072
255 7.350424 7.543097
256 7.280742 7.475080
257 7.457409 7.657053
258 7.362589 7.562601
259 7.488655 7.691595
260 7.397386 7.600527
261 7.496311 7.701572
262 7.415863 7.621325
263 7.496871 7.704075
264 7.428099 7.635545
265 7.495639 7.704599
266 7.437547 7.646801
267 7.494288 7.704897
268 7.445469 7.656424
269 7.493305 7.705494
270 7.452385 7.664969
271 7.492793 7.706513
272 7.458553 7.672712


[[16]]
[[16]]$mean
Time Series:
Start = 192 
End = 211 
Frequency = 1 
 [1] 1.580294 1.556735 1.540904 1.526367 1.512048 1.497764 1.483487 1.469211 1.454934 1.440658 1.426382 1.412106 1.397830 1.383554 1.369278 1.355002
[17] 1.340726 1.326450 1.312174 1.297898

[[16]]$lower
Time Series:
Start = 192 
End = 211 
Frequency = 1 
          80%       95%
192 1.2781168 1.1181537
193 1.2269430 1.0523619
194 1.1952168 1.0122211
195 1.1668251 0.9764950
196 1.1393757 0.9420951
197 1.1124437 0.9084673
198 1.0859255 0.8754692
199 1.0597751 0.8430330
200 1.0339597 0.8111091
201 1.0084524 0.7796563
202 0.9832296 0.7486387
203 0.9582708 0.7180249
204 0.9335578 0.6877868
205 0.9090742 0.6578997
206 0.8848057 0.6283415
207 0.8607391 0.5990922
208 0.8368628 0.5701338
209 0.8131660 0.5414500
210 0.7896390 0.5130258
211 0.7662729 0.4848478

[[16]]$upper
Time Series:
Start = 192 
End = 211 
Frequency = 1 
         80%      95%
192 1.882472 2.042435
193 1.886526 2.061107
194 1.886591 2.069586
195 1.885909 2.076240
196 1.884719 2.082000
197 1.883085 2.087061
198 1.881048 2.091505
199 1.878646 2.095388
200 1.875909 2.098760
201 1.872864 2.101660
202 1.869535 2.104126
203 1.865942 2.106188
204 1.862103 2.107874
205 1.858034 2.109208
206 1.853750 2.110214
207 1.849265 2.110912
208 1.844589 2.111318
209 1.839734 2.111450
210 1.834708 2.111322
211 1.829522 2.110947


[[17]]
[[17]]$mean
Time Series:
Start = 224 
End = 243 
Frequency = 1 
 [1] 234.1930 234.6295 234.3313 234.2910 234.2454 234.2271 234.2165 234.2113 234.2085 234.2071 234.2063 234.2060 234.2058 234.2057 234.2056 234.2056
[17] 234.2056 234.2056 234.2055 234.2055

[[17]]$lower
Time Series:
Start = 224 
End = 243 
Frequency = 1 
         80%      95%
224 185.9335 160.3865
225 183.2916 156.1149
226 181.1893 153.0577
227 180.3707 151.8270
228 179.8300 151.0243
229 179.4525 150.4566
230 179.1418 149.9870
231 178.8660 149.5680
232 178.6078 149.1745
233 178.3590 148.7947
234 178.1155 148.4228
235 177.8752 148.0555
236 177.6371 147.6915
237 177.4005 147.3297
238 177.1653 146.9699
239 176.9311 146.6119
240 176.6980 146.2553
241 176.4659 145.9003
242 176.2347 145.5467
243 176.0044 145.1946

[[17]]$upper
Time Series:
Start = 224 
End = 243 
Frequency = 1 
         80%      95%
224 282.4525 307.9996
225 285.9674 313.1441
226 287.4732 315.6049
227 288.2114 316.7551
228 288.6607 317.4664
229 289.0017 317.9976
230 289.2912 318.4460
231 289.5565 318.8545
232 289.8092 319.2425
233 290.0552 319.6194
234 290.2972 319.9899
235 290.5367 320.3564
236 290.7744 320.7201
237 291.0108 321.0816
238 291.2459 321.4413
239 291.4800 321.7993
240 291.7131 322.1558
241 291.9452 322.5108
242 292.1764 322.8644
243 292.4067 323.2165


[[18]]
[[18]]$mean
Time Series:
Start = 224 
End = 243 
Frequency = 1 
 [1] 2.293675 2.335371 2.311209 2.311209 2.311209 2.311209 2.311209 2.311209 2.311209 2.311209 2.311209 2.311209 2.311209 2.311209 2.311209 2.311209
[17] 2.311209 2.311209 2.311209 2.311209

[[18]]$lower
Time Series:
Start = 224 
End = 243 
Frequency = 1 
          80%        95%
224 0.6232249 -0.2610580
225 0.6164538 -0.2934861
226 0.5473775 -0.3863385
227 0.5353904 -0.4046712
228 0.5234836 -0.4228810
229 0.5116557 -0.4409703
230 0.4999050 -0.4589414
231 0.4882300 -0.4767968
232 0.4766293 -0.4945385
233 0.4651015 -0.5121687
234 0.4536453 -0.5296895
235 0.4422593 -0.5471029
236 0.4309422 -0.5644109
237 0.4196928 -0.5816153
238 0.4085100 -0.5987180
239 0.3973924 -0.6157208
240 0.3863391 -0.6326254
241 0.3753490 -0.6494334
242 0.3644208 -0.6661466
243 0.3535537 -0.6827664

[[18]]$upper
Time Series:
Start = 224 
End = 243 
Frequency = 1 
         80%      95%
224 3.964126 4.848409
225 4.054289 4.964228
226 4.075040 5.008756
227 4.087027 5.027089
228 4.098934 5.045299
229 4.110762 5.063388
230 4.122513 5.081359
231 4.134188 5.099215
232 4.145789 5.116956
233 4.157316 5.134587
234 4.168773 5.152107
235 4.180159 5.169521
236 4.191476 5.186829
237 4.202725 5.204033
238 4.213908 5.221136
239 4.225025 5.238139
240 4.236079 5.255043
241 4.247069 5.271851
242 4.257997 5.288564
243 4.268864 5.305184


[[19]]
[[19]]$mean
Time Series:
Start = 197 
End = 216 
Frequency = 1 
 [1] 0.3410592 0.2961105 0.3348482 0.3014633 0.3302351 0.3054389 0.3268088 0.3083918 0.3242640 0.3105850 0.3223738 0.3122139 0.3209699 0.3134238
[15] 0.3199272 0.3143225 0.3191528 0.3149899 0.3185776 0.3154857

[[19]]$lower
Time Series:
Start = 197 
End = 216 
Frequency = 1 
           80%        95%
197 -0.5323444 -0.9946963
198 -0.6089233 -1.0880192
199 -0.5779045 -1.0610865
200 -0.6376093 -1.1347243
201 -0.6181518 -1.1201976
202 -0.6656152 -1.1796602
203 -0.6547316 -1.1743278
204 -0.6932252 -1.2234491
205 -0.6886550 -1.2248619
206 -0.7205224 -1.2663577
207 -0.7205704 -1.2726717
208 -0.7475077 -1.3084904
209 -0.7509101 -1.3183291
210 -0.7741537 -1.3498825
211 -0.7799733 -1.3622254
212 -0.8004298 -1.3905440
213 -0.8079736 -1.4046382
214 -0.8263109 -1.4304791
215 -0.8350684 -1.4457717
216 -0.8517805 -1.4696939

[[19]]$upper
Time Series:
Start = 197 
End = 216 
Frequency = 1 
         80%      95%
197 1.214463 1.676815
198 1.201144 1.680240
199 1.247601 1.730783
200 1.240536 1.737651
201 1.278622 1.780668
202 1.276493 1.790538
203 1.308349 1.827945
204 1.310009 1.840233
205 1.337183 1.873390
206 1.341692 1.887528
207 1.365318 1.917419
208 1.371936 1.932918
209 1.392850 1.960269
210 1.401001 1.976730
211 1.419828 2.002080
212 1.429075 2.019189
213 1.446279 2.042944
214 1.456291 2.060459
215 1.472224 2.082927
216 1.482752 2.100665

#ALL DATA INTERPOLATED

interpolationfunction=function(k){
  data3=lakedata%>%
    filter(characteristic_name== k)%>%
    group_by(year)%>%
    summarise(mean=mean(result_value))
  allyears=seq(1982,2022,1)
  dfallyears=data.frame(year=allyears)
  merged_data <- merge(dfallyears, data3, by = "year", all = TRUE)
  finaldata=merged_data%>%
    select(mean)
  na_p_data2 <- na.interp(finaldata)
  added_data_P=as.data.frame(na_p_data2)
  autoplot(na_p_data2, series="Interpolated") +
    autolayer(as.ts(finaldata), series="Original") +
    geom_point()+
    scale_colour_manual(values=c(`Interpolated`="red",`Original`="Blue"))+
    labs(title = k) +geom_vline(xintercept = c(22, 26), linetype = "dashed", color = "red")
}

interplatedlist=list()
for( i in 1:(length(characteristic.names)-1)){
   plot.data= interpolationfunction(characteristic.names[i])
 interplatedlist[[i]] <- plot.data
}
interplatedlist

[[1]]

[[2]]

[[3]]

[[4]]

[[5]]

Warning: Removed 11 rows containing missing values (`geom_line()`).

[[6]]

[[7]]

[[8]]

[[9]]

[[10]]

[[11]]

[[12]]

[[13]]

[[14]]

[[15]]

[[16]]

[[17]]

[[18]]

[[19]]

Warning: Removed 7 rows containing missing values (`geom_line()`).